#ifndef LATTICE_H
#define LATTICE_H
#include <vector>
#include "point.h"
#include <algorithm>
#include <stdio.h>

namespace{
    #define sqr(x) ((x)*(x))
    constexpr double PI=3.1415926535897932384626433832795;
    const char*  const Symbol_A ="A";
    const char*  const Symbol_B ="B";
    const char*  const Symbol_C ="C";
    const char*  const Symbol_D ="D";
    const char*  const Symbol_E ="E";
    const char*  const Symbol_F ="F";
    const char*  const Symbol_G ="Γ";
    const char*  const Symbol_H ="H";
    const char*  const Symbol_K ="K";
    const char*  const Symbol_L ="L";
    const char*  const Symbol_M ="M";
    const char*  const Symbol_N ="N";
    const char*  const Symbol_P ="P";
    const char*  const Symbol_Q ="Q";
    const char*  const Symbol_R ="R";
    const char*  const Symbol_S ="S";
    const char*  const Symbol_T ="T";
    const char*  const Symbol_U ="U";
    const char*  const Symbol_W ="W";
    const char*  const Symbol_X ="X";
    const char*  const Symbol_Y ="Y";
    const char*  const Symbol_Z ="Z";
    const char*  const Symbol_0 ="Ω";
    const char*  const Symbol_A1="A_{1}";
    const char*  const Symbol_B1="B_{1}";
    const char*  const Symbol_C1="C_{1}";
    const char*  const Symbol_D1="D_{1}";
    const char*  const Symbol_H1="H_{1}";
    const char*  const Symbol_H2="H_{2}";
    const char*  const Symbol_L1="L_{1}";
    const char*  const Symbol_L2="L_{2}";
    const char*  const Symbol_M1="M_{1}";
    const char*  const Symbol_M2="M_{2}";
    const char*  const Symbol_N1="N_{1}";
    const char*  const Symbol_P1="P_{1}";
    const char*  const Symbol_P2="P_{2}";
    const char*  const Symbol_Q1="Q_{1}";
    const char*  const Symbol_S1="Σ";
    const char*  const Symbol_S2="Σ_{1}";
    const char*  const Symbol_X1="X_{1}";
    const char*  const Symbol_Y1="Y_{1}";
    const char*  const Symbol_Y2="Y_{2}";
    const char*  const Symbol_Y3="Y_{3}";
    const char*  const Symbol_Z1="Z_{1}";
//typedef unsigned int uint;
}

template <typename real>
class Latt
{
private:
    point<real> data[3];
    std::vector<point<real>> HSP;
    std::vector<const char*> symbol;
    int    type;
public:
    static real  ERROR;
    Latt(){type=0;}
    int  GetType(){
        if(! type) FinSym();
        return type;
    }
    void operator=(const Latt &orig){
        memcpy(data,orig.data,sizeof(real)*9);
    }
    void operator=(real num){
        for(int i=0;i<3;i++){
            for(int j=0;j<3;j++){
                data[i][j]=num;
            }
        }
    }
    point<real> &operator[](const int &i){
        return data[i];
    }
    point<real> operator[](const int &i)const {
        return data[i];
    }
    real &operator()(const int &i,const int &j){
        return data[i][j];
    }
    real  operator()(const int &i,const int &j) const {
        return data[i][j];
    }
    void  operator*=(real A){
        data[0]*=A;data[1]*=A;data[2]*=A;
    }
    Latt  operator*(const Latt &N){
        Latt<real> T;
        for(int i=0;i<3;i++){
            for(int k=0;k<3;k++){
                for(int j=0;j<3;j++){
                    T.data[i][j]+=data[i][k]*N.data[k][j];
                }
            }
        }
        return T;
    }
    Latt  &operator*=(const Latt &N){
        Latt<real> T;
        for(int i=0;i<3;i++){
            for(int k=0;k<3;k++){
                for(int j=0;j<3;j++){
                    T.data[i][j]+=data[i][k]*N.data[k][j];
                }
            }
        }
        *this=T;
        return *this;
    }    
    void  swapcol(const int a,const int b){
        for(int i=0;i<3;i++){
            std::swap(data[i][a],data[i][b]);
        }
    }
    void  swaprow(const int a,const int b){
        for(int i=0;i<3;i++){
            std::swap(data[a][i],data[b][i]);
        }
    }
    void  Cart(){        
        Latt<real>RLV=GetRLV();
        for(int i=0;i<18;i++){
            HSP[i]*=RLV;
        }
    }
    real  CosAgl(const int a,const int b){
        return data[a] * data[b] / (data[a].norm()*data[b].norm() );
    }
    Latt  GetRLV(){        
        Latt<real> RLV=inv(*this);
        RLV*=(2*PI);
        std::swap(RLV[0][1],RLV[1][0]);
        std::swap(RLV[0][2],RLV[2][0]);
        std::swap(RLV[1][2],RLV[2][1]);
        return RLV;
    }
    void  HSPSymbol();
    void  FindHSP();
    void  ShowHSP(int n=6){
        Latt<real>RLV_1=inv(GetRLV());
        for(int i=0;i<static_cast<int>(HSP.size());i++){
            HSP[i].show(n);
            if(symbol[i]){
                std::cout<<' '<<symbol[i];
            }
            (HSP[i]*RLV_1).show(n);
            std::cout<<std::endl;
        }
    }
    Latt<real> FinSym();
    std::vector<int> StdLatt(){
        std::vector<int> swap_order;
        if(data[0].norm()>data[1].norm()){
            swapcol(0,1);
            swaprow(0,1);
            swap_order.push_back(0);
            swap_order.push_back(1);
        }
        if(data[1].norm()>data[2].norm()){
            swapcol(1,2);
            swaprow(1,2);
            swap_order.push_back(1);
            swap_order.push_back(2);
        }
        if(data[0].norm()>data[1].norm()){
            swapcol(0,1);
            swaprow(0,1);
            swap_order.push_back(0);
            swap_order.push_back(1);
        }
        while(1){
            real x;
            while(1){
                x=data[0].norm();
                x=round(data[0]*data[1]/(x*x));
                data[1]=data[1]-x*data[0];
                if(data[0].norm()<=data[1].norm()){
                    break;
                }else{
                    swapcol(0,1);
                    swaprow(0,1);
                    swap_order.push_back(0);
                    swap_order.push_back(1);
                }
            }
            x=data[1].norm();
            x=round(data[1]*data[2]/(x*x));
            data[2]=data[2]-x*data[1];
            if(data[1].norm()<=data[2].norm()){
                break;
            }else{
                swapcol(1,2);
                swaprow(1,2);
                swap_order.push_back(1);
                swap_order.push_back(2);
            }
        }
        /*
        while(swap_order.size()){
            int a=swap_order.back();
            swap_order.pop_back();
            int b=swap_order.back();
            swap_order.pop_back();
            swapcol(a,b);
            swaprow(a,b);
        }
        */
       return swap_order;
    }
    const char* HighSymPnt(const point<real> &p);
};

template <typename real>
void  point<real>::operator*=(const Latt<real>  &M)
{
    point<real> V=*this;
    P[0]=V[0]*M[0][0]+V[1]*M[1][0]+V[2]*M[2][0];
    P[1]=V[0]*M[0][1]+V[1]*M[1][1]+V[2]*M[2][1];
    P[2]=V[0]*M[0][2]+V[1]*M[1][2]+V[2]*M[2][2];
}

template <typename real>
point<real> operator*(const point<real> &V,const Latt<real> &M)
{
    point<real> T;
    T[0]=V[0]*M[0][0]+V[1]*M[1][0]+V[2]*M[2][0];
    T[1]=V[0]*M[0][1]+V[1]*M[1][1]+V[2]*M[2][1];
    T[2]=V[0]*M[0][2]+V[1]*M[1][2]+V[2]*M[2][2];
    return T;
}

template <typename real>
point<real> operator*(const Latt<real> &M,const point<real> &V)
{
    point<real> T;
    T[0]=M[0][0]*V[0]+M[0][1]*V[1]+M[0][2]*V[2];
    T[1]=M[1][0]*V[0]+M[1][1]*V[1]+M[1][2]*V[2];
    T[2]=M[2][0]*V[0]+M[2][1]*V[1]+M[2][2]*V[2];
    return T;
}

template <typename real>
Latt<real> inv(const Latt<real> &M){
    
    Latt<real> I;
    real const bM=fabs( M[0] * cross(M[1],M[2]) );
    //printf("%lf\n",bM);
    I(0,0)=( M(1,1) * M(2,2) - M(1,2) * M(2,1) ) / bM;
    I(1,0)=( M(1,2) * M(2,0) - M(1,0) * M(2,2) ) / bM;
    I(2,0)=( M(1,0) * M(2,1) - M(1,1) * M(2,0) ) / bM;
    I(0,1)=( M(0,2) * M(2,1) - M(0,1) * M(2,2) ) / bM;
    I(1,1)=( M(0,0) * M(2,2) - M(0,2) * M(2,0) ) / bM;
    I(2,1)=( M(0,1) * M(2,0) - M(0,0) * M(2,1) ) / bM;
    I(0,2)=( M(0,1) * M(1,2) - M(0,2) * M(1,1) ) / bM;
    I(1,2)=( M(0,2) * M(1,0) - M(0,0) * M(1,2) ) / bM;
    I(2,2)=( M(0,0) * M(1,1) - M(0,1) * M(1,0) ) / bM;
    return I;

}

template <typename real>
Latt<real> RotMat(const point<real> &axis,const real &angle)
{
    Latt<real> M;
    M[0][0]= cos(angle);
    M[0][1]=-sin(angle);
    M[0][2]= 0;
    M[1][0]= sin(angle);
    M[1][1]= cos(angle);
    M[1][2]= 0;
    M[2][0]= 0;
    M[2][1]= 0;
    M[2][2]= 1;
    if(sqrt(axis[0]*axis[0]+axis[1]*axis[1])/fabs(axis[2])<1E-10){
        return M;
    }
    real theta=asin(axis[2]/axis.norm());
    real phi  =atan(axis[1]/axis[0]);
    if(axis[0]<0) phi  +=PI;
    if(axis[1]<0) theta+=PI;
    Latt<real> P,T;
    P[0][0]= cos(-phi);
    P[0][1]=-sin(-phi);
    P[0][2]= 0;
    P[1][0]= sin(-phi);
    P[1][1]= cos(-phi);
    P[1][2]= 0;
    P[2][0]= 0;
    P[2][1]= 0;
    P[2][2]= 1;

    T[0][0]= cos(theta);
    T[0][1]= 0;
    T[0][2]= sin(theta);
    T[1][0]= 0;
    T[1][1]= 1;
    T[1][2]= 0;
    T[2][0]=-sin(theta);
    T[2][1]= 0;
    T[2][2]= cos(theta);

    M=inv(P)*inv(T)*M*T*P;
    return M;
}

template <typename real>
real Latt<real>::ERROR=1E-5;



template <typename real>
Latt<real> Latt<real>::FinSym()
{
    Latt<real> DLV=*this;
    std::vector<int>swap;
    //DLV.StdLatt();
    for(int i=0;i<3;i++){ //冒泡排序从小到大排
        for(int j=i+1;j<3;j++){
            if(DLV[i].norm()>DLV[j].norm()){
                DLV.swapcol(i,j);
                DLV.swaprow(i,j);
                swap.push_back(i);
                swap.push_back(j);                
            }
        }
    }
    real A=DLV[0].norm();
    real B=DLV[1].norm();
    real C=DLV[2].norm();
    real alpha=fabs(DLV.CosAgl(1,2));
    real beta =fabs(DLV.CosAgl(0,2));
    real gamma=fabs(DLV.CosAgl(0,1));
    if( fabs(alpha-beta) < ERROR && fabs(alpha-gamma)< ERROR){//α=β=γ
        if(alpha<ERROR){//α=β=γ=90°
            if(fabs(A-B)<ERROR && fabs(A-C)<ERROR ){// A=B=C立方晶系 Cubic
                type= 1;
                goto END;
            }else if(fabs(A-B)<ERROR ){
                type=4;
                goto END;
            }else if(fabs(B-C)<ERROR ){
                DLV.swapcol(0,2);
                DLV.swaprow(0,2);
                swap.push_back(0);
                swap.push_back(2);
                type=4;
                goto END;
            }else{
                type=7; //Orthorhombic 正交
                goto END;
            }
        }else if(fabs(alpha- 0.5)<ERROR && fabs(A-B)<ERROR && fabs(A-C)<ERROR ){//α=60° 面心立方 Face-centered cubic
            type= 2;
            goto END;
        }else if(fabs(alpha-1/3.)<ERROR && fabs(A-B)<ERROR && fabs(A-C)<ERROR ){//α=109°28' 体心立方 Body-centered cubic
            type= 3;
            goto END;
        }else if(fabs(A-B)<ERROR && fabs(A-C)<ERROR ){ //三方晶系 Rhombohedral α=β=γ  A=B=C
            alpha=DLV.CosAgl(1,2);
            beta =DLV.CosAgl(0,2);
            gamma=DLV.CosAgl(0,1);

            real angle;
            if(fabs(alpha-beta)<ERROR){
                angle=gamma;
            }else if(fabs(alpha-gamma)<ERROR){
                angle=beta;
            }else{
                angle=alpha;
            }
            if(angle>0){//Rhombohedral
                type= 14;
                goto END;
            }else{
                type= 15;
                goto END;
            }
        }
    }else if(fabs(A-B)<ERROR && fabs(A-C)<ERROR){// only A=B=C
        if(fabs(alpha-gamma) < ERROR){//α=γ
            DLV.swapcol(0,2);  //α=β
            DLV.swaprow(0,2);
            swap.push_back(0);
            swap.push_back(2);
        }else if(fabs(beta-gamma) < ERROR){//β=γ
            DLV.swapcol(1,2);  //α=β
            DLV.swaprow(1,2);
            swap.push_back(1);
            swap.push_back(2);
        }//let α=β
        alpha=fabs(DLV.CosAgl(1,2));
        beta =fabs(DLV.CosAgl(0,2));
        gamma=fabs(DLV.CosAgl(0,1));
        if(fabs(alpha-beta) < ERROR){
            if(alpha<ERROR){//α=β=90° 
                if(fabs(gamma-0.5) <ERROR){//γ=120°
                    type= 13;
                    goto END;
                }else{
                    type= 12;
                    goto END;
                }
            }else{//体心四方 Body-centered tetragonal            
                alpha=acos(alpha);
                A=2*A*sin(alpha/2);//得到惯用胞中晶格常数A的数值 A/2/L=sin(α/2), L为元胞晶格长度
                C=fabs(cross(DLV[0],DLV[1])*DLV[2])/A/A*2;
                if(A>C){
                    type= 5;
                    goto END;
                }else{
                    type= 6;
                    goto END;
                }

            }
        }else{//体心正交 α>β>γ Body-centered orthorhombic
            type= 11;
            goto END;
        }
    }else if(alpha+beta<ERROR || alpha+gamma <ERROR || beta+gamma <ERROR){ //某一个向量垂直于另外两个向量组成的平面
        if(alpha+gamma <ERROR){ //γ=90° A⊥B α=90° B⊥C
            DLV.swaprow(1,2);
            DLV.swapcol(1,2);//C⊥AB
            swap.push_back(1);
            swap.push_back(2);
        }else if(beta+gamma <ERROR){ //γ=90° A⊥B β=90° A⊥C 
            DLV.swaprow(0,2);
            DLV.swapcol(0,2);//C⊥AB
            swap.push_back(0);
            swap.push_back(2);
        }//在二维材料中，通常c轴作为真空层方向，因此调整为C⊥AB，而不是传统的A⊥BC
        if(DLV[0].norm()>DLV[1].norm()){ //let A<B
            DLV.swaprow(0,1);
            DLV.swapcol(0,1);
            swap.push_back(0);
            swap.push_back(1);
        }
        A=DLV[0].norm();
        B=DLV[1].norm();
        C=DLV[2].norm();
        alpha=fabs(DLV.CosAgl(1,2));
        beta =fabs(DLV.CosAgl(0,2));
        gamma=fabs(DLV.CosAgl(0,1));

        if(fabs(A-B)<ERROR && fabs(gamma-1/2.)<ERROR){//γ=120° or 60° Hexagonal
            type=13;
            goto END;
        }else if(fabs(A-B)<ERROR){//C-centered orthorhombic
            type= 12;
            goto END;

        }else{
            type= 16; //Monoclinic 单斜
            goto END;
        }
    }
    {
        real const sqr_c=2*(B*B-A*A+C*C);
        real const sqr_b=2*(A*A-B*B+C*C);
        real const sqr_a=2*(A*A+B*B-C*C);
        real const Va=sqrt(sqr_a*sqr_b*sqr_c);    
        real const Vb=fabs(cross(DLV[0],DLV[1])*DLV[2]);
        if(fabs(Va-4*Vb)<ERROR){ //Face-centered orhtorhombic 面心正交
            if(fabs(1/sqr_a-1/sqr_b-1/sqr_c) <ERROR){
                type= 9;
                goto END;
            }else if(1/sqr_a > 1/sqr_b+1/sqr_c){
                type= 8;
                goto END;
            }else{
                type= 10;
                goto END;
            }
        }
    }

    if(fabs(B-C)<ERROR){
        DLV.swaprow(0,2);
        DLV.swapcol(0,2);
        swap.push_back(0);
        swap.push_back(2);        
    }//let A=B
    A=DLV[0].norm();
    B=DLV[1].norm();
    C=DLV[2].norm();
    alpha=DLV.CosAgl(1,2);
    beta =DLV.CosAgl(0,2);
    gamma=DLV.CosAgl(0,1);
    if(fabs(A-B)<ERROR){ //C-centered monoclinic 底心单斜
        Latt<real> RLV=DLV.GetRLV();
        real kgamma=RLV.CosAgl(0,1);
        if(fabs(kgamma)<ERROR){
            type= 18;
            goto END;
        }else if(kgamma<0){
            type= 17;
            goto END;
        }
        real AdotB=data[0]*data[1];
        real ConLattB=sqrt(2*(A+AdotB));
        real ConLattA=sqrt(2*(A-AdotB));
        real TanAlpha=data[2][2]/data[2][1];
        real Alpha=atan(TanAlpha);
        //real ConLattC=data[2][2]/sin(Alpha);
        if(fabs(ConLattB*cos(Alpha)/C+sqr(ConLattB)*sqr(sin(Alpha))/sqr(ConLattA)-1)<ERROR){
            type= 20;
            goto END;
        }else if(ConLattB*cos(Alpha)/C+sqr(ConLattB)*sqr(sin(Alpha))/sqr(ConLattA)<1){
            type= 19;
            goto END;
        }else{
            type= 21;
            goto END;
        }
    }

    //Triclinic 三斜
    {
        Latt<real> RLV=DLV.GetRLV();
        real kalpha=RLV.CosAgl(1,2);
        real kbeta =RLV.CosAgl(0,2);
        real kgamma=RLV.CosAgl(0,1);
        if(fabs(kgamma)<ERROR && kalpha<0 && kbeta <0){
            type= 23;
            goto END;
        }else if(fabs(kgamma)<ERROR && kalpha>0 && kbeta >0){
            type= 25;
            goto END;
        }else if(kalpha<0 && kbeta <0 && kgamma<0){
            type= 22;
            goto END;
        }else if(kalpha>0 && kbeta >0 && kgamma>0){
            type= 24;
            goto END;
        }
    }
END:
    return DLV;
}

template <typename real>
void  Latt<real>::HSPSymbol()
{
    if(!HSP.size()) FindHSP();
    Latt<real> DLV=FinSym();
    Latt<real> RLV;
    symbol.back()=Symbol_G;
    switch(type){
        case 1:{
            point<real> M(1/2., 1/2.,   0 );
            point<real> R(1/2., 1/2., 1/2.);
            point<real> X(  0 , 1/2.,   0 );
            DLV=0;
            DLV[0][0]=1;
            DLV[1][1]=1;
            DLV[2][2]=1;
            DLV*=data[0].norm();
            RLV =DLV.GetRLV();
            M*=RLV;
            R*=RLV;
            X*=RLV;
            real const LM=M.norm();
            real const LR=R.norm();
            real const LX=X.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LM)<ERROR){
                    symbol[i]=Symbol_M;
                    continue;
                }
                if(fabs(HSP[i].norm()-LR)<ERROR){
                    symbol[i]=Symbol_R;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                }
            }

            break;
        }
        case 2:{
            point<real> K(3/8., 3/8., 3/4.);
            point<real> L(1/2., 1/2., 1/2.);
            point<real> U(5/8., 1/4., 5/8.);
            point<real> W(1/2., 1/4., 3/4.);
            point<real> X(1/2.,   0 , 1/2.);
            DLV=1;
            DLV[0][0]=0;
            DLV[1][1]=0;
            DLV[2][2]=0;
            DLV*=data[0].norm()/sqrt(2.);
            RLV=DLV.GetRLV();
            K*=RLV;
            L*=RLV;
            U*=RLV;
            W*=RLV;
            X*=RLV;
            real const LK=K.norm();
            real const LL=L.norm();
            real const LW=W.norm();
            real const LX=X.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LK)<ERROR){
                    symbol[i]=Symbol_K;
                    continue;
                }
                if(fabs(HSP[i].norm()-LL)<ERROR){
                    symbol[i]=Symbol_L;
                    continue;
                }
                if(fabs(HSP[i].norm()-LW)<ERROR){
                    symbol[i]=Symbol_W;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                    continue;
                }
            }
            real const UX=(U-X).norm();
            int startX=0,endX=0;
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i] == Symbol_X){
                    startX=i;
                    break;
                }
            }
            for(int i=startX;i<static_cast<int>(HSP.size());i++){
                if(symbol[i] != Symbol_X){
                    endX=i;
                    break;
                }
            }
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i] != Symbol_K) continue;
                for(int j=startX;j<endX;j++){
                    if(fabs((HSP[i]-HSP[j]).norm()-UX)<ERROR){
                        symbol[i]=Symbol_U;
                        break;
                    }
                }
            }

            break;
        }
        case 3:{
            point<real> N(   0 ,    0., 1/2.);
            point<real> H( 1/2., -1/2., 1/2.);
            point<real> P( 1/4.,  1/4., 1/4.);
            DLV=1;
            DLV[0][0]*=-1;
            DLV[1][1]*=-1;
            DLV[2][2]*=-1;
            DLV*=data[0].norm()/sqrt(3.);
            RLV=DLV.GetRLV();
            N*=RLV;
            H*=RLV;
            P*=RLV;
            real const LN=N.norm();
            real const LH=H.norm();
            real const LP=P.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LN)<ERROR){
                    symbol[i]=Symbol_N;
                    continue;
                }
                if(fabs(HSP[i].norm()-LH)<ERROR){
                    symbol[i]=Symbol_H;
                    continue;
                }
                if(fabs(HSP[i].norm()-LP)<ERROR){
                    symbol[i]=Symbol_P;
                }
            }
            break;
        }
        case 4:{
            real const LLA=DLV[0].norm();
            real const LLC=DLV[2].norm();
            point<real> A(1/2., 1/2., 1/2.);
            point<real> M(1/2., 1/2.,   0 );
            point<real> R(  0 , 1/2., 1/2.);
            point<real> X(  0 , 1/2.,   0 );
            point<real> Z(  0 ,   0 , 1/2.);
            
            DLV=0;
            DLV[0][0]=LLA;
            DLV[1][1]=LLA;
            DLV[2][2]=LLC;
            RLV=DLV.GetRLV();
            A*=RLV;
            M*=RLV;
            R*=RLV;
            X*=RLV;
            Z*=RLV;
            real const LA=A.norm();
            real const LM=M.norm();
            real const LR=R.norm();
            real const LX=X.norm();
            real const LZ=Z.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LA)<ERROR){
                    symbol[i]=Symbol_A;
                    continue;
                }
                if(fabs(HSP[i].norm()-LM)<ERROR){
                    symbol[i]=Symbol_M;
                    continue;
                }
                if(fabs(HSP[i].norm()-LR)<ERROR){
                    symbol[i]=Symbol_R;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                    continue;
                }
                if(fabs(HSP[i].norm()-LZ)<ERROR){
                    symbol[i]=Symbol_Z;
                    continue;
                }
            }
            break;
        }
        case 5:{
            real const alpha=acos(fabs(DLV.CosAgl(1,2)));//得到alpha的绝对值
            real const A=2*DLV[0].norm()*sin(alpha/2);//得到惯用胞中晶格常数A的数值 A/2/L=sin(α/2), L为元胞晶格长度
            real const C=cross(DLV[0],DLV[1])*DLV[2]/A/A*2;
            real const eta =(1 + sqr(C) / sqr(A) )/4;

            DLV[0][0]=-A/2;
            DLV[0][1]= A/2;
            DLV[0][2]= C/2;
            DLV[1][0]= A/2;
            DLV[1][1]=-A/2;
            DLV[1][2]= C/2;
            DLV[2][0]= A/2;
            DLV[2][1]= A/2;
            DLV[2][2]=-C/2;
            RLV=DLV.GetRLV();

            point<real> M (-1/2.,  1/2., 1/2.);
            point<real> N (   0 ,  1/2.,   0 );
            point<real> P ( 1/4.,  1/4., 1/4.);
            point<real> X (   0 ,    0 , 1/2.);
            point<real> Z ( eta ,  eta ,-eta );
            point<real> Z1(-eta ,1-eta , eta );
            M *=RLV;
            N *=RLV;
            P *=RLV;
            X *=RLV;
            Z *=RLV;
            Z1*=RLV;
            real const LM =M.norm();
            real const LN =N.norm();
            real const LP =P.norm();
            real const LX =X.norm();
            real const LZ =Z.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LM)<ERROR){
                    symbol[i]=Symbol_M;
                    continue;
                }
                if(fabs(HSP[i].norm()-LN)<ERROR){
                    symbol[i]=Symbol_N;
                    continue;
                }
                if(fabs(HSP[i].norm()-LP)<ERROR){
                    symbol[i]=Symbol_P;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                    continue;
                }
                if(fabs(HSP[i].norm()-LZ)<ERROR){
                    symbol[i]=Symbol_Z;
                    continue;
                }
            }
            int startM=0,endM=0;
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]==Symbol_M){
                    startM=i;
                    while(symbol[i]==Symbol_M) i++;
                    endM=i;
                    break;
                }
            }
            real const MZ1=(Z1-M).norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_Z) continue;
                for(int j=startM;j<endM;j++){
                    if(fabs((HSP[i]-HSP[j]).norm()-MZ1)<ERROR){
                        symbol[i]=Symbol_Z1;
                        break;
                    }
                }
            }
            break;
        }
        case 6:{
            real const alpha=acos(fabs(DLV.CosAgl(1,2)));//得到alpha的绝对值
            real const A=2*DLV[0].norm()*sin(alpha/2);//得到惯用胞中晶格常数A的数值 A/2/L=sin(α/2), L为元胞晶格长度
            real const C=cross(DLV[0],DLV[1])*DLV[2]/A/A*2;
            real const eta =(1+ sqr(A)/sqr(C))/4;
            real const zeta=sqr(A)/(2*sqr(C));

            DLV[0][0]=-A/2;
            DLV[0][1]= A/2;
            DLV[0][2]= C/2;
            DLV[1][0]= A/2;
            DLV[1][1]=-A/2;
            DLV[1][2]= C/2;
            DLV[2][0]= A/2;
            DLV[2][1]= A/2;
            DLV[2][2]=-C/2;
            RLV=DLV.GetRLV();

            point<real> N (    0 ,  1/2.,    0 );
            point<real> P (  1/4.,  1/4.,  1/4.);
            point<real> S ( -eta ,  eta ,  eta );//Σ
            point<real> S1(  eta ,1-eta , -eta );//Σ1
            point<real> X (    0 ,    0 ,  1/2.);
            point<real> Y (-zeta , zeta ,  1/2.);
            point<real> Y1(  1/2.,  1/2., -zeta);
            point<real> Z (  1/2.,  1/2., -1/2.);
            N *=RLV;
            P *=RLV;
            S *=RLV;
            S1*=RLV;
            X *=RLV;
            Y *=RLV;
            Y1*=RLV;
            Z *=RLV;
            
            real const LN=N.norm();
            real const LP=P.norm();
            real const LS=S.norm();
            real const LX=X.norm();
            real const LY=Y.norm();
            real const LZ=Z.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LN)<ERROR){
                    symbol[i]=Symbol_N;
                    continue;
                }
                if(fabs(HSP[i].norm()-LP)<ERROR){
                    symbol[i]=Symbol_P;
                    continue;
                }
                if(fabs(HSP[i].norm()-LS)<ERROR){
                    symbol[i]=Symbol_S1;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                    continue;
                }
                if(fabs(HSP[i].norm()-LY)<ERROR){
                    symbol[i]=Symbol_Y;
                    continue;
                }
                if(fabs(HSP[i].norm()-LZ)<ERROR){
                    symbol[i]=Symbol_Z;
                    continue;
                }
            }
            int startZ=0,endZ=0;
            for(int i=0;i<static_cast<int>(symbol.size());i++){
                if(symbol[i]==Symbol_Z){
                    startZ=i;
                    while(symbol[i] == Symbol_Z) i++;
                    endZ=i;
                    break;
                }
            }
            real const ZS1=(Z-S1).norm();
            real const ZY1=(Z-Y1).norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_S1) continue;
                for(int j=startZ;j<endZ;j++){
                    if(fabs((HSP[i]-HSP[j]).norm()-ZS1)<ERROR){
                        symbol[i]=Symbol_S2;
                        break;
                    }
                }
            }
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_Y) continue;
                for(int j=startZ;j<endZ;j++){
                    if(fabs((HSP[i]-HSP[j]).norm()-ZY1)<ERROR){
                        symbol[i]=Symbol_Y1;
                        break;
                    }
                }
            }
            break;
        }
        case 7:{
            real const A=DLV[0].norm();
            real const B=DLV[1].norm();
            real const C=DLV[2].norm();
            DLV=0;
            DLV[0][0]=A;
            DLV[1][1]=B;
            DLV[2][2]=C;
            RLV=DLV.GetRLV();

            point<real> R(1/2.,1/2.,1/2.);
            point<real> S(1/2.,1/2.,  0 );
            point<real> T(  0 ,1/2.,1/2.);
            point<real> U(1/2.,  0 ,1/2.);
            point<real> X(1/2.,  0 ,  0 );
            point<real> Y(  0 ,1/2.,  0 );
            point<real> Z(  0 ,  0 ,1/2.);
            R*=RLV;
            S*=RLV;
            T*=RLV;
            U*=RLV;
            X*=RLV;
            Y*=RLV;
            Z*=RLV;

            real const LR=R.norm();
            real const LS=S.norm();
            real const LT=T.norm();
            real const LU=U.norm();
            real const LX=X.norm();
            real const LY=Y.norm();
            real const LZ=Z.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LR)<ERROR){
                    symbol[i]=Symbol_R;
                    continue;
                }
                if(fabs(HSP[i].norm()-LS)<ERROR){
                    symbol[i]=Symbol_S;
                    continue;
                }
                if(fabs(HSP[i].norm()-LT)<ERROR){
                    symbol[i]=Symbol_T;
                    continue;
                }
                if(fabs(HSP[i].norm()-LU)<ERROR){
                    symbol[i]=Symbol_U;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                    continue;
                }
                if(fabs(HSP[i].norm()-LY)<ERROR){
                    symbol[i]=Symbol_Y;
                    continue;
                }
                if(fabs(HSP[i].norm()-LZ)<ERROR){
                    symbol[i]=Symbol_Z;
                    continue;
                }
            }
            break;
        }
        case 8:{
            real LA=DLV[0].norm();
            real LB=DLV[1].norm();
            real LC=DLV[2].norm();
            real const sqr_c=2*(LB*LB-LA*LA+LC*LC);
            real const sqr_b=2*(LA*LA-LB*LB+LC*LC);
            real const sqr_a=2*(LA*LA+LB*LB-LC*LC);
            real const zeta=(1+sqr_a/sqr_b-sqr_a/sqr_c)/4;
            real const  eta=(1+sqr_a/sqr_b+sqr_a/sqr_c)/4;
            LA=sqrt(sqr_a);
            LB=sqrt(sqr_b);
            LC=sqrt(sqr_c);
            DLV[0][0]=   0;
            DLV[0][1]=LB/2;
            DLV[0][2]=LC/2;
            DLV[1][0]=LA/2;
            DLV[1][1]=   0;
            DLV[1][2]=LC/2;
            DLV[2][0]=LA/2;
            DLV[2][1]=LB/2;
            DLV[2][2]=   0;
            RLV=DLV.GetRLV();
            point<real> A (1/2., 1/2. + zeta ,   zeta );
            point<real> A1(1/2., 1/2. - zeta , 1-zeta );
            point<real> L (1/2.,         1/2.,    1/2.);
            point<real> T (  1 ,         1/2.,    1/2.);
            point<real> X (  0 ,         eta ,    eta );
            point<real> X1(  1 ,       1-eta ,  1-eta );
            point<real> Y (1/2.,           0 ,    1/2.);
            point<real> Z (1/2.,         1/2.,      0 );
            A *=RLV;
            A1*=RLV;
            L *=RLV;
            T *=RLV;
            X *=RLV;
            X1*=RLV;
            Y *=RLV;
            Z *=RLV;
                       LA=A.norm(); 
            real const LL=L.norm(); 
            real const LT=T.norm(); 
            real const LX=X.norm(); 
            real const LY=Y.norm(); 
            real const LZ=Z.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LA)<ERROR){
                    symbol[i]=Symbol_A;
                    continue;
                }
                if(fabs(HSP[i].norm()-LL)<ERROR){
                    symbol[i]=Symbol_L;
                    continue;
                }
                if(fabs(HSP[i].norm()-LT)<ERROR){
                    symbol[i]=Symbol_T;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                    continue;
                }
                if(fabs(HSP[i].norm()-LY)<ERROR){
                    symbol[i]=Symbol_Y;
                    continue;
                }
                if(fabs(HSP[i].norm()-LZ)<ERROR){
                    symbol[i]=Symbol_Z;
                    continue;
                }
            }
            
            int startZ=0,endZ=0;
            for(int i=0;i<static_cast<int>(symbol.size());i++){
                if(symbol[i]==Symbol_Z){
                    startZ=i;
                    while(symbol[i]==Symbol_Z) i++;
                    endZ=i;
                    break;
                }
            }

            real const ZX1=(Z-X1).norm();
            real const ZA1=(Z-A1).norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_X) continue;
                for(int j=startZ;j<endZ;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-ZX1 ) < ERROR ){
                        symbol[i]=Symbol_X1;
                        break;
                    }
                }
            }
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_A) continue;
                for(int j=startZ;j<endZ;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-ZA1 ) < ERROR ){
                        symbol[i]=Symbol_A1;
                        break;
                    }
                }
            }

            break;
        }
        case 9:{
            real LA=DLV[0].norm();
            real LB=DLV[1].norm();
            real LC=DLV[2].norm();
            real const sqr_c=2*(LB*LB-LA*LA+LC*LC);
            real const sqr_b=2*(LA*LA-LB*LB+LC*LC);
            real const sqr_a=2*(LA*LA+LB*LB-LC*LC);
            real const zeta=(1+sqr_a/sqr_b-sqr_a/sqr_c)/4;
            real const  eta=(1+sqr_a/sqr_b+sqr_a/sqr_c)/4;
            LA=sqrt(sqr_a);
            LB=sqrt(sqr_b);
            LC=sqrt(sqr_c);
            DLV[0][0]=   0;
            DLV[0][1]=LB/2;
            DLV[0][2]=LC/2;
            DLV[1][0]=LA/2;
            DLV[1][1]=   0;
            DLV[1][2]=LC/2;
            DLV[2][0]=LA/2;
            DLV[2][1]=LB/2;
            DLV[2][2]=   0;
            RLV=DLV.GetRLV();
            point<real> A (1/2., 1/2. + zeta ,   zeta );
            point<real> A1(1/2., 1/2. - zeta , 1-zeta );
            point<real> L (1/2.,         1/2.,    1/2.);
            point<real> T (  1 ,         1/2.,    1/2.);
            point<real> X (  0 ,         eta ,    eta );
            point<real> Y (1/2.,           0 ,    1/2.);
            point<real> Z (1/2.,         1/2.,      0 );
            A *=RLV;
            A1*=RLV;
            L *=RLV;
            T *=RLV;
            X *=RLV;
            Y *=RLV;
            Z *=RLV;
                       LA=A.norm(); 
            real const LL=L.norm(); 
            real const LT=T.norm(); 
            real const LX=X.norm(); 
            real const LY=Y.norm(); 
            real const LZ=Z.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LA)<ERROR){
                    symbol[i]=Symbol_A;
                    continue;
                }
                if(fabs(HSP[i].norm()-LL)<ERROR){
                    symbol[i]=Symbol_L;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                    continue;
                }
                if(fabs(HSP[i].norm()-LT)<ERROR){
                    symbol[i]=Symbol_T;
                    continue;
                }
                if(fabs(HSP[i].norm()-LY)<ERROR){
                    symbol[i]=Symbol_Y;
                    continue;
                }
                if(fabs(HSP[i].norm()-LZ)<ERROR){
                    symbol[i]=Symbol_Z;
                    continue;
                }
            }
            
            int startZ=0,endZ=0;
            for(int i=0;i<static_cast<int>(symbol.size());i++){
                if(symbol[i]==Symbol_Z){
                    startZ=i;
                    while(symbol[i]==Symbol_Z) i++;
                    endZ=i;
                    break;
                }
            }

            real const ZT =(Z-T ).norm();
            real const ZA1=(Z-A1).norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_X) continue;
                for(int j=startZ;j<endZ;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-ZT ) < ERROR ){
                        symbol[i]=Symbol_T;
                        break;
                    }
                }
            }
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_A) continue;
                for(int j=startZ;j<endZ;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-ZA1 ) < ERROR ){
                        symbol[i]=Symbol_A1;
                        break;
                    }
                }
            }

            break;
        }
        case 10:{
            real LA=DLV[0].norm();
            real LB=DLV[1].norm();
            real LC=DLV[2].norm();
            real const sqr_c=2*(LB*LB-LA*LA+LC*LC);
            real const sqr_b=2*(LA*LA-LB*LB+LC*LC);
            real const sqr_a=2*(LA*LA+LB*LB-LC*LC);
            real const eta  =(1+sqr_a/sqr_b-sqr_a/sqr_c)/4;
            real const phi  =(1+sqr_c/sqr_b-sqr_c/sqr_a)/4;
            real const delta=(1+sqr_b/sqr_a-sqr_b/sqr_c)/4;
            LA=sqrt(sqr_a);
            LB=sqrt(sqr_b);
            LC=sqrt(sqr_c);
            DLV[0][0]=   0;
            DLV[0][1]=LB/2;
            DLV[0][2]=LC/2;
            DLV[1][0]=LA/2;
            DLV[1][1]=   0;
            DLV[1][2]=LC/2;
            DLV[2][0]=LA/2;
            DLV[2][1]=LB/2;
            DLV[2][2]=   0;
            RLV=DLV.GetRLV();

            point<real> C (         1/2., 1/2. - eta ,     1-eta );
            point<real> C1(         1/2., 1/2. + eta ,       eta );
            point<real> D (1/2. - delta ,        1/2., 1 - delta );
            point<real> D1(1/2. + delta ,        1/2.,     delta );
            point<real> L (         1/2.,        1/2.,       1/2.);
            point<real> H (       1-phi , 1/2. - phi ,       1/2.);
            point<real> H1(         phi , 1/2. + phi ,       1/2.);
            point<real> X (           0 ,        1/2.,       1/2.);
            point<real> Y (         1/2.,          0 ,       1/2.);
            point<real> Z (         1/2.,        1/2.,         0 );
            C *=RLV;
            C1*=RLV;
            D *=RLV;
            D1*=RLV;
            L *=RLV;
            H *=RLV;
            H1*=RLV;
            X *=RLV;
            Y *=RLV;
            Z *=RLV;

                       LC =C.norm();
            real const LD =D.norm();
            real const LL =L.norm();
            real const LH =H.norm();
            real const LX =X.norm();
            real const LY =Y.norm();
            real const LZ =Z.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LC)<ERROR){
                    symbol[i]=Symbol_C;
                    continue;
                }
                if(fabs(HSP[i].norm()-LD)<ERROR){
                    symbol[i]=Symbol_D;
                    continue;
                }
                if(fabs(HSP[i].norm()-LL)<ERROR){
                    symbol[i]=Symbol_L;
                    continue;
                }
                if(fabs(HSP[i].norm()-LH)<ERROR){
                    symbol[i]=Symbol_H;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                    continue;
                }
                if(fabs(HSP[i].norm()-LY)<ERROR){
                    symbol[i]=Symbol_Y;
                    continue;
                }
                if(fabs(HSP[i].norm()-LZ)<ERROR){
                    symbol[i]=Symbol_Z;
                    continue;
                }
            }

            int startZ=0,endZ=0;
            for(int i=0;i<static_cast<int>(symbol.size());i++){
                if(symbol[i]==Symbol_Z){
                    startZ=i;
                    while(symbol[i]==Symbol_Z) i++;
                    endZ=i;
                    break;
                }
            }
            int startX=0,endX=0;
            for(int i=0;i<static_cast<int>(symbol.size());i++){
                if(symbol[i]==Symbol_X){
                    startX=i;
                    while(symbol[i]==Symbol_X) i++;
                    endX=i;
                    break;
                }
            }

            real const ZC1=(Z-C1).norm();
            real const ZD1=(Z-D1).norm();
            real const XH1=(X-H1).norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_C) continue;
                for(int j=startZ;j<endZ;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-ZC1 ) < ERROR ){
                        symbol[i]=Symbol_C1;
                        break;
                    }
                }
            }
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_D) continue;
                for(int j=startZ;j<endZ;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-ZD1 ) < ERROR ){
                        symbol[i]=Symbol_D1;
                        break;
                    }
                }
            }
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_H) continue;
                for(int j=startX;j<endX;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-XH1 ) < ERROR ){
                        symbol[i]=Symbol_H1;
                        break;
                    }
                }
            }
            break;
        }
        case 11:{//α>β>γ cos(γ)>cos(β)>cos(α)
            
            real A=DLV[0].norm();
            real B=DLV[1].norm();
            real C=DLV[2].norm();
            real alpha=fabs(DLV.CosAgl(1,2));
            real beta =fabs(DLV.CosAgl(0,2));
            real gamma=fabs(DLV.CosAgl(0,1));

            if(gamma<beta ) std::swap(gamma,beta);
            if(gamma<alpha) std::swap(gamma,alpha);
            if(alpha>beta ) std::swap(alpha,beta);

            gamma=acos(gamma);
            A=2*A* sin(gamma/2);  //A/2/L=sin(angle/2)
            beta =acos(beta);
            B=2*B* sin(beta/2);
            alpha=acos(-alpha);
            C=2*C* sin( alpha/2);
            DLV[0][0]=-A/2;
            DLV[0][1]= B/2;
            DLV[0][2]= C/2;
            DLV[1][0]= A/2;
            DLV[1][1]=-B/2;
            DLV[1][2]= C/2;
            DLV[2][0]= A/2;
            DLV[2][1]= B/2;
            DLV[2][2]=-C/2;
            RLV=DLV.GetRLV();
            

            real const AdotB=DLV[0]*DLV[1];//c^2/4-a^2/4-b^2/4
            real const AdotC=DLV[0]*DLV[2];//b^2/4-a^2/4-c^2/4
            real const BdotC=DLV[1]*DLV[2];//a^2/4-b^2/4-c^2/4
            real const sqra =2*fabs(AdotB+AdotC);
            real const sqrb =2*fabs(AdotB+BdotC);
            real const sqrc =2*fabs(AdotC+BdotC);
            real const zeta =(1+sqra/sqrc)/4;
            real const eta  =(1+sqrb/sqrc)/4;
            real const delta=(sqrb-sqra)/(4*sqrc);
            real const miu  =(sqra+sqrb)/(4*sqrc);

            point<real> L (         -miu ,          miu , 1/2. - delta );
            point<real> L1(          miu ,         -miu , 1/2. + delta );
            point<real> L2( 1/2. - delta , 1/2. + delta ,         -miu );
            point<real> R (            0 ,          1/2.,            0 );
            point<real> S (          1/2.,            0 ,            0 );
            point<real> T (            0 ,            0 ,          1/2.);
            point<real> W (          1/4.,          1/4.,          1/4.);
            point<real> X (        -zeta ,         zeta ,         zeta );
            point<real> X1(         zeta ,     1 - zeta ,        -zeta );
            point<real> Y (          eta ,         -eta ,          eta );
            point<real> Y1(      1 - eta ,          eta ,         -eta );
            point<real> Z (          0.5 ,          0.5 ,         -0.5 );
            L *=RLV;
            L1*=RLV;
            L2*=RLV;
            R *=RLV;
            S *=RLV;
            T *=RLV;
            W *=RLV;
            X *=RLV;
            X1*=RLV;
            Y *=RLV;
            Y1*=RLV;
            Z *=RLV;
            
            real const LL=L.norm();
            real const LR=R.norm();
            real const LS=S.norm();
            real const LT=T.norm();
            real const LW=W.norm();
            real const LX=X.norm();
            real const LY=Y.norm();
            real const LZ=Z.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LL)<ERROR){
                    symbol[i]=Symbol_L;
                    continue;
                }
                if(fabs(HSP[i].norm()-LR)<ERROR){
                    symbol[i]=Symbol_R;
                    continue;
                }
                if(fabs(HSP[i].norm()-LS)<ERROR){
                    symbol[i]=Symbol_S;
                    continue;
                }
                if(fabs(HSP[i].norm()-LT)<ERROR){
                    symbol[i]=Symbol_T;
                    continue;
                }
                if(fabs(HSP[i].norm()-LW)<ERROR){
                    symbol[i]=Symbol_W;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                    continue;
                }
                if(fabs(HSP[i].norm()-LY)<ERROR){
                    symbol[i]=Symbol_Y;
                    continue;
                }
                if(fabs(HSP[i].norm()-LZ)<ERROR){
                    symbol[i]=Symbol_Z;
                    continue;
                }
            }

            int startZ=0,endZ=0;
            for(int i=0;i<static_cast<int>(symbol.size());i++){
                if(symbol[i]==Symbol_Z){
                    startZ=i;
                    while(symbol[i]==Symbol_Z) i++;
                    endZ=i;
                    break;
                }
            }

            int startY=0,endY=0;
            for(int i=0;i<static_cast<int>(symbol.size());i++){
                if(symbol[i]==Symbol_Y){
                    startY=i;
                    while(symbol[i]==Symbol_Y) i++;
                    endY=i;
                    break;
                }
            }
            real const ZX1=(Z-X1).norm();
            real const ZY1=(Z-Y1).norm();
            real const ZL2=(Z-L2).norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_X) continue;
                for(int j=startZ;j<endZ;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-ZX1 ) < ERROR ){
                        symbol[i]=Symbol_X1;
                        break;
                    }
                }
            }
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_Y) continue;
                for(int j=startZ;j<endZ;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-ZY1 ) < ERROR ){
                        symbol[i]=Symbol_Y1;
                        break;
                    }
                }
            }
            
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_L) continue;
                for(int j=startZ;j<endZ;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-ZL2 ) < ERROR ){
                        symbol[i]=Symbol_L2;
                        break;
                    }
                }
            }
            
            real const YL1=(Y-L1).norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_L) continue;
                for(int j=startY;j<endY;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-YL1 ) < ERROR ){
                        symbol[i]=Symbol_L1;
                        break;
                    }
                }
            }
            if(fabs(sqr(C) - sqr(B) -sqr(A) )<ERROR){
                int startS=0,endS=0;
                for(int i=0;i<static_cast<int>(symbol.size());i++){
                    if(symbol[i]==Symbol_S){
                        startS=i;
                        while(symbol[i]==Symbol_S) i++;
                        endS=i;
                        break;
                    }
                }
            
                real const SW=(S-W).norm();
                for(int i=0;i<static_cast<int>(HSP.size());i++){
                    if(symbol[i]!=Symbol_L) continue;
                    for(int j=startS;j<endS;j++){
                        if( fabs( (HSP[i]-HSP[j]).norm()-SW ) < ERROR ){
                            symbol[i]=Symbol_W;
                            break;
                        }
                    }
                }
            }

            break;
        }
        case 12:{
            real LA=DLV[0].norm();
            real LB=DLV[1].norm();
            real LC=DLV[2].norm();
            real const AdotB=fabs(DLV[0]*DLV[1]);
            real const sqrb =(LA*LA+AdotB)*2;
            real const sqra =(LA*LA-AdotB)*2;
            real const zeta =(1+sqra/sqrb)/4;
            LB=sqrt(sqrb);
            LA=sqrt(sqra);
            DLV=0;
            DLV[0][0]= LA/2;
            DLV[0][1]=-LB/2;
            DLV[1][0]= LA/2;
            DLV[1][1]= LB/2;
            DLV[2][2]= LC;
            RLV=DLV.GetRLV();

            point<real> A ( zeta ,     zeta , 1/2.);
            point<real> A1( -zeta, 1 - zeta , 1/2.);
            point<real> R (    0 ,      1/2., 1/2.);
            point<real> S (    0 ,      1/2.,   0 );
            point<real> T ( -1/2.,      1/2., 1/2.);
            point<real> X ( zeta ,     zeta ,   0 );
            point<real> X1(-zeta ,   1-zeta ,   0 );
            point<real> Y ( -1/2.,      1/2.,   0 );
            point<real> Z (    0 ,        0 , 1/2.);
            A *=RLV;
            A1*=RLV;
            R *=RLV;
            S *=RLV;
            T *=RLV;
            X *=RLV;
            X1*=RLV;
            Y *=RLV;
            Z *=RLV;

                       LA=A.norm();
            real const LR=R.norm();
            real const LS=S.norm();
            real const LT=T.norm();
            real const LX=X.norm();
            real const LY=Y.norm();
            real const LZ=Z.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LA)<ERROR){
                    symbol[i]=Symbol_A;
                    continue;
                }
                if(fabs(HSP[i].norm()-LR)<ERROR){
                    symbol[i]=Symbol_R;
                    continue;
                }
                if(fabs(HSP[i].norm()-LS)<ERROR){
                    symbol[i]=Symbol_S;
                    continue;
                }
                if(fabs(HSP[i].norm()-LT)<ERROR){
                    symbol[i]=Symbol_T;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                    continue;
                }
                if(fabs(HSP[i].norm()-LY)<ERROR){
                    symbol[i]=Symbol_Y;
                    continue;
                }
                if(fabs(HSP[i].norm()-LZ)<ERROR){
                    symbol[i]=Symbol_Z;
                    continue;
                }
            }

            int startY=0,endY=0;
            for(int i=0;i<static_cast<int>(symbol.size());i++){
                if(symbol[i]==Symbol_Y){
                    startY=i;
                    break;
                }
            }
            for(int i=startY;i<static_cast<int>(symbol.size());i++){
                if(symbol[i] != Symbol_Y){
                    endY=i;
                    break;
                }
            }

            real const YA1=(Y-A1).norm();
            real const YX1=(Y-X1).norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_A) continue;
                for(int j=startY;j<endY;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-YA1 ) < ERROR ){
                        symbol[i]=Symbol_A1;
                        break;
                    }
                }
            }
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_X) continue;
                for(int j=startY;j<endY;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-YX1 ) < ERROR ){
                        symbol[i]=Symbol_X1;
                        break;
                    }
                }
            }

            break;
        }
        case 13:{
            real LA=DLV[0].norm();
            real LC=DLV[2].norm();
            DLV=0;
            DLV[0][0]= LA;
            DLV[1][0]=-LA/2;
            DLV[1][1]=(LA/2)*sqrt(3.);
            DLV[2][2]= LC;
            RLV=DLV.GetRLV();
            point<real> A(  0.,  0.,1/2.);
            point<real> H(1/3.,1/3.,1/2.);
            point<real> K(1/3.,1/3.,  0.);
            point<real> L(1/2.,  0.,1/2.);
            point<real> M(1/2.,  0.,  0.);
            A*=RLV;
            H*=RLV;
            K*=RLV;
            L*=RLV;
            M*=RLV;
            
                       LA=A.norm();
            real const LH=H.norm();
            real const LK=K.norm();
            real const LL=L.norm();
            real const LM=M.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LA)<ERROR){
                    symbol[i]=Symbol_A;
                    continue;
                }
                if(fabs(HSP[i].norm()-LH)<ERROR){
                    symbol[i]=Symbol_H;
                    continue;
                }
                if(fabs(HSP[i].norm()-LK)<ERROR){
                    symbol[i]=Symbol_K;
                    continue;
                }
                if(fabs(HSP[i].norm()-LL)<ERROR){
                    symbol[i]=Symbol_L;
                    continue;
                }
                if(fabs(HSP[i].norm()-LM)<ERROR){
                    symbol[i]=Symbol_M;
                    continue;
                }
            }

            break;
        }
        case 14:{
            real A=DLV[0].norm();
            real alpha=DLV.CosAgl(1,2);
            real beta =DLV.CosAgl(0,2);
            real gamma=DLV.CosAgl(0,1);

            real angle;
            if(fabs(alpha-beta)<ERROR){
                angle=gamma;
            }else if(fabs(alpha-gamma)<ERROR){
                angle=beta;
            }else{
                angle=alpha;
            }
            real const eta  =(1+4*angle)/(2+4*angle);
            real const niu  =3/4.-eta/2;
            angle=acos(angle);
            DLV=0;
            DLV[0][0]= A*cos(angle/2.);
            DLV[0][1]=-A*sin(angle/2.);
            DLV[1][0]= A*cos(angle/2.);
            DLV[1][1]= A*sin(angle/2.);
            DLV[2][0]= A*cos(angle)/cos(angle/2);
            DLV[2][2]= A*sqrt( 1 - sqr(cos(angle))/sqr(cos(angle/2)) );
            RLV=DLV.GetRLV();

            point<real> B (    eta ,     1/2., 1 - eta );
            point<real> B1(    1/2., 1 - eta , eta - 1 );
            point<real> F (    1/2.,     1/2.,       0 );
            point<real> L (    1/2.,       0 ,       0 );
            point<real> L1(      0 ,       0 ,    -1/2.);
            point<real> P (    eta ,     niu ,     niu );
            point<real> P1(1 - niu , 1 - niu , 1 - eta );
            point<real> P2(    niu ,     niu , eta - 1 );
            point<real> Q (1 - niu ,     niu ,       0 );
            point<real> X (    niu ,       0 ,    -niu );
            point<real> Z (    1/2.,     1/2.,     1/2.);
            B *=RLV;
            B1*=RLV;
            F *=RLV;
            L *=RLV;
            L1*=RLV;
            P *=RLV;
            P1*=RLV;
            P2*=RLV;
            Q *=RLV;
            X *=RLV;
            Z *=RLV;

            real const LB=B.norm();
            real const LF=F.norm();
            real const LL=L.norm();
            real const LP=P.norm();
            real const LQ=Q.norm();
            real const LX=X.norm();
            real const LZ=Z.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LB)<ERROR){
                    symbol[i]=Symbol_B;
                    continue;
                }
                if(fabs(HSP[i].norm()-LF)<ERROR){
                    symbol[i]=Symbol_F;
                    continue;
                }
                if(fabs(HSP[i].norm()-LL)<ERROR){
                    symbol[i]=Symbol_L;
                    continue;
                }
                if(fabs(HSP[i].norm()-LP)<ERROR){
                    symbol[i]=Symbol_P;
                    continue;
                }
                if(fabs(HSP[i].norm()-LQ)<ERROR){
                    symbol[i]=Symbol_Q;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                    continue;
                }
                if(fabs(HSP[i].norm()-LZ)<ERROR){
                    symbol[i]=Symbol_Z;
                    continue;
                }
            }
            int startF=0,endF=0;
            for(int i=0;i<static_cast<int>(symbol.size());i++){
                if(symbol[i]==Symbol_F){
                    startF=i;
                    while(symbol[i]==Symbol_F) i++;
                    endF=i;
                    break;
                }
            }
            real const FQ=(F-Q).norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_Q) continue;
                for(int j=startF;j<endF;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-FQ ) > ERROR ){
                        symbol[i]=Symbol_X;
                        break;
                    }
                }
            }

            int startL=0,endL=0;
            for(int i=0;i<static_cast<int>(symbol.size());i++){
                if(symbol[i]==Symbol_L){
                    startL=i;
                    while(symbol[i]==Symbol_L) i++;
                    endL=i;
                    break;
                }
            }
            real const LP1=(L-P1).norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_P) continue;
                for(int j=startL;j<endL;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-LP1 ) < ERROR ){
                        symbol[i]=Symbol_P1;
                        break;
                    }
                }
            }
            int startZ=0,endZ=0;
            for(int i=0;i<static_cast<int>(symbol.size());i++){
                if(symbol[i]==Symbol_Z){
                    startZ=i;
                    while(symbol[i]==Symbol_Z) i++;
                    endZ=i;
                    break;
                }
            }
            real const ZB1=(Z-B1).norm();
            real const ZP2=(Z-P2).norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_B) continue;
                for(int j=startZ;j<endZ;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-ZB1 ) < ERROR ){
                        symbol[i]=Symbol_B1;
                        break;
                    }
                }
            }
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_P1) continue;
                for(int j=startZ;j<endZ;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-ZP2 ) < ERROR ){
                        symbol[i]=Symbol_P2;
                        break;
                    }
                }
            }
            break;
        }
        case 15:{            
            real A=DLV[0].norm();
            real alpha=DLV.CosAgl(1,2);
            real beta =DLV.CosAgl(0,2);
            real gamma=DLV.CosAgl(0,1);

            real angle;
            if(fabs(alpha-beta)<ERROR){
                angle=gamma;
            }else if(fabs(alpha-gamma)<ERROR){
                angle=beta;
            }else{
                angle=alpha;
            }
            angle=acos(angle);

            DLV=0;
            DLV[0][0]= A*cos(angle/2.);
            DLV[0][1]=-A*sin(angle/2.);
            DLV[1][0]= A*cos(angle/2.);
            DLV[1][1]= A*sin(angle/2.);
            DLV[2][0]= A*cos(angle)/cos(angle/2);
            DLV[2][2]= A*sqrt(1-sqr(cos(angle))/sqr(cos(angle/2)) );
            RLV=DLV.GetRLV();
            real const eta  =1/( 2*sqr( tan(angle/2) ) );
            real const niu  =3/4.-eta/2;

            point<real> F (    1/2.,    -1/2.,       0 );
            point<real> L (    1/2.,       0 ,       0 );
            point<real> P (1 - niu ,    -niu , 1 - niu );
            point<real> P1(    niu , niu - 1 , niu - 1 );
            point<real> Q (    eta ,     eta ,     eta );
            point<real> Q1(1 - eta ,    -eta ,    -eta );
            point<real> Z (    1/2.,    -1/2.,     1/2.);

            F *=RLV;
            L *=RLV;
            P *=RLV;
            P1*=RLV;
            Q *=RLV;
            Q1*=RLV;
            Z *=RLV;
        
            real const LF=F.norm();
            real const LL=L.norm();
            real const LP=P.norm();
            real const LQ=Q.norm();
            real const LZ=Z.norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LF)<ERROR){
                    symbol[i]=Symbol_F;
                    continue;
                }
                if(fabs(HSP[i].norm()-LL)<ERROR){
                    symbol[i]=Symbol_L;
                    continue;
                }
                if(fabs(HSP[i].norm()-LP)<ERROR){
                    symbol[i]=Symbol_P;
                    continue;
                }
                if(fabs(HSP[i].norm()-LQ)<ERROR){
                    symbol[i]=Symbol_Q;
                    continue;
                }
                if(fabs(HSP[i].norm()-LZ)<ERROR){
                    symbol[i]=Symbol_Z;
                    continue;
                }
            }

            int startF=0,endF=0;
            for(int i=0;i<static_cast<int>(symbol.size());i++){
                if(symbol[i]==Symbol_F){
                    startF=i;
                    while(symbol[i]==Symbol_F) i++;
                    endF=i;
                    break;
                }
            }

            real const FQ1=(F-Q1).norm();
            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(symbol[i]!=Symbol_Q) continue;
                for(int j=startF;j<endF;j++){
                    if( fabs( (HSP[i]-HSP[j]).norm()-FQ1 ) < ERROR ){
                        symbol[i]=Symbol_Q1;
                        break;
                    }
                }
            }
            //缺少判别P与P1点的方法
            break;
        }
        case 16:{// c⊥ab a<b
            real LA=DLV[0].norm();
            real LB=DLV[1].norm();
            real LC=DLV[2].norm();
            real const gamma=acos(fabs(DLV.CosAgl(0,1)));
            real const eta=(1-LA*cos(gamma)/LB)/( 2 * sqr( sin(gamma) ) );
            real const miu=1/2.-eta*LB*cos(gamma)/LA;
            DLV=0;
            DLV[0][0]=LA;
            DLV[2][2]=LC;
            DLV[1][0]=LB*cos(gamma);
            DLV[1][1]=LB*sin(gamma);
            RLV=DLV.GetRLV();
            
            point<real> A (   1/2.,      0 , 1/2.);
            point<real> C (   1/2.,    1/2.,   0 );
            point<real> D (     0 ,    1/2., 1/2.);
            point<real> D1(     0 ,  - 1/2., 1/2.);
            point<real> E (   1/2.,    1/2., 1/2.);
            point<real> H (    eta,1 - miu ,   0 );
            point<real> H1(1 - eta,    miu ,   0 );
            point<real> H2(    eta,  - miu ,   0 );
            point<real> M (    eta,1 - miu , 1/2.);
            point<real> M1(1 - eta,    miu , 1/2.);
            point<real> M2(    eta,  - miu , 1/2.);
            point<real> X (   1/2.,      0 ,   0 );
            point<real> Y (     0 ,    1/2.,   0 );
            point<real> Y1(     0 ,  - 1/2.,   0 );
            point<real> Z (     0 ,      0 , 1/2.);
            A *=RLV;
            C *=RLV;
            D *=RLV;
            D1*=RLV;
            E *=RLV;
            H *=RLV;
            H1*=RLV;
            H2*=RLV;
            M *=RLV;
            M1*=RLV;
            M2*=RLV;
            X *=RLV;
            Y *=RLV;
            Y1*=RLV;
            Z *=RLV;
            
                       LA=A.norm();
                       LC=C.norm();
            real const LD=D.norm();
            real const LE=E.norm();
            real const LH=H.norm();
            real const LM=M.norm();
            real const LX=X.norm();
            real const LY=Y.norm();
            real const LZ=Z.norm();

            for(int i=0;i<static_cast<int>(HSP.size());i++){
                if(fabs(HSP[i].norm()-LA)<ERROR){
                    symbol[i]=Symbol_A;
                    continue;
                }
                if(fabs(HSP[i].norm()-LC)<ERROR){
                    symbol[i]=Symbol_C;
                    continue;
                }
                if(fabs(HSP[i].norm()-LD)<ERROR){
                    symbol[i]=Symbol_D;
                    continue;
                }
                if(fabs(HSP[i].norm()-LE)<ERROR){
                    symbol[i]=Symbol_E;
                    continue;
                }
                if(fabs(HSP[i].norm()-LH)<ERROR){
                    symbol[i]=Symbol_H;
                    continue;
                }
                if(fabs(HSP[i].norm()-LM)<ERROR){
                    symbol[i]=Symbol_M;
                    continue;
                }
                if(fabs(HSP[i].norm()-LX)<ERROR){
                    symbol[i]=Symbol_X;
                    continue;
                }
                if(fabs(HSP[i].norm()-LY)<ERROR){
                    symbol[i]=Symbol_Y;
                    continue;
                }
                if(fabs(HSP[i].norm()-LZ)<ERROR){
                    symbol[i]=Symbol_Z;
                    continue;
                }
            }

            break;
        }
        case 17:
            break;
        case 18:
            break;
        case 19:
            break;
        case 20:
            break;
        case 21:
            break;
        case 22:
            break;
        case 23:
            break;
        case 24:
            break;
        case 25:
            break;
    }
    for(int i=0;i<static_cast<int>(symbol.size());i++){
        if(!symbol[i]){
            symbol[i]=symbol.back();
            HSP[i]=HSP.back();
            symbol.pop_back();
            HSP.pop_back();
            i--;
        }
    }
}

template <typename real>
const char* Latt<real>::HighSymPnt(const point<real> &p)
{
    for(int i=0;i<static_cast<int>(HSP.size());i++){
        if( (p-HSP[i]).norm()<ERROR ){
            return symbol[i];
        }
    }
    return NULL;
}

template <typename real>
void Latt<real>::FindHSP()
{
    Latt<real> RLV=GetRLV();
    point<real> Vec[26];
    int counter=0;
    for(int i=-1;i<2;i++){
        for(int j=-1;j<2;j++){
            for(int k=-1;k<2;k++){
                if(i==0 && j==0 && k==0) continue;
                Vec[counter]=(RLV[0]*i+RLV[1]*j+RLV[2]*k)/2.;
                counter++;
            }
        }
    }


    std::vector<real> planed;
    for(int i=0;i<26;i++){
        HSP.push_back(Vec[i]);
        planed.push_back(Vec[i]*Vec[i]);
    }

    Latt<real>  Matrix;//记录由三个非共面向量组成的矩阵
    point<real> vecd;

    //printf("%20.16lf%20.16lf%20.16lf %lf\n",Vec[26][0],Vec[26][1],Vec[26][2],planed[26]);
    //过点vec[i]且垂直于向量vec[i]的平面方程为 vec[i][0]*x+vec[i][1]*y+vec[i][2]*z-vec[i]*vec[i]=0
    for(int i=0;i<counter;i++){
        Matrix[0]=Vec[i];
        vecd[0]=planed[i];
        for(int j=i+1;j<counter;j++){
            Matrix[1]=Vec[j];
            vecd[1]=planed[j];
            for(int k=j+1;k<counter;k++){
                if(fabs(Matrix[0] * cross(Matrix[1],Vec[k]) ) <ERROR) continue;
                Matrix[2]=Vec[k];
                vecd[2]=planed[k];
                HSP.push_back(inv(Matrix)*vecd);
            }
        }
    }


    for(int i=0;i<static_cast<int>(HSP.size());i++){
        for(int j=0;j<26;j++){
            real E=Vec[j]*HSP[i]-planed[j];
            if(-planed[j] * E < -ERROR){
                HSP[i]=HSP.back();
                HSP.pop_back();
                i--;
                break;
            }
        }
    }

    //删除重复的点
    for(int i=0;i<static_cast<int>(HSP.size());i++){
        for(int j=i+1;j<static_cast<int>(HSP.size());j++){
            if( (HSP[i]-HSP[j]).norm() < ERROR){
                HSP[j]=HSP.back();
                HSP.pop_back();
                j--;
            }
        }
    }

    std::vector<point<real>*> cross_point; //记录三面交点
    for(int i=0;i<static_cast<int>(HSP.size());i++){
        counter=0;
        for(int j=0;j<26;j++){
            real E=Vec[j]*HSP[i]-planed[j];
            if( fabs(E)  < ERROR) counter++;
            if(counter==3){
                cross_point.push_back(&HSP[i]);
                break;
            }
        }
    }
    //添加任意两个交点的中点，并判断其是否在棱上
    for(int i=0;i<static_cast<int>(cross_point.size());i++){
        for(int j=i+1;j<static_cast<int>(cross_point.size());j++){            
            point<real> middle=(*cross_point[i]+*cross_point[j])/2.;
            counter=0;
            for(int k=0;k<26;k++){
                real E=Vec[k]*middle-planed[k];
                if( fabs(E)  < ERROR) counter++;
            }
            if(counter==2) HSP.push_back(middle);
        }
    }
    //删除重复的点
    for(int i=0;i<static_cast<int>(HSP.size());i++){
        for(int j=i+1;j<static_cast<int>(HSP.size());j++){
            if( (HSP[i]-HSP[j]).norm() < ERROR){
                HSP[j]=HSP.back();
                HSP.pop_back();
                j--;
            }
        }
    }
    std::sort(HSP.begin(), HSP.end(), [](const point<real> &a, const point<real> &b) {
        return a.norm() > b.norm();
    });
    /*
    for(int i=0;i<static_cast<int>(HSP.size());i++){
        HSP[i].show();
    }
    std::cout<<std::endl;
    */
    point<real> gamma;
    HSP.push_back(gamma);
    symbol.resize(static_cast<int>(HSP.size()));
    symbol.assign(static_cast<int>(symbol.size()),0);

}

#endif